% Version of calibrateaNM_byHand.m for model with free entry.
% See Model with Free Entry Appendix

function [w,sM,sN, varrho,Timplied] = calibrateaNM_byHand_FE(sD, sX, phi,fc,fx,param)

eta = param.eta;
g = param.gamma;
eps = param.eps;
sg = param.sg;
beta = param.beta;
qD = param.qD;
theta = param.theta;
vartheta = param.vartheta;
zeta = param.zeta;
pM = param.pM;
AN = param.AN;
FXM = param.FXM; FX = param.FX; FM = param.FM;

% Constants

k1 =  qD^(g*(sg-1)) * (1-g)^((1-g)*(sg-1)) * g^(g*(sg-1)) * g * (sg-1)^sg/sg^sg * eta;  
k2 = qD^(g*(sg-1)) * (sg-1)^(sg-1)/sg^sg * (1-g)^((1-g)*(sg-1))* g^(g*(sg-1));
k3 = k2 * (theta-1)/theta;
k5 = (sg / (sg-1))^(1-sg) * (1-g)^((1-g)*(sg-1)) * (g*beta^(eps/(eps-1))*qD)^(g*(sg-1)); 
k6 = mean(k1*(sD.^(-1)-1).^(1/(eps-1)/eta) .* fc + k2 .* ((sD<1).*(sX==0)) * FM + k3 * (sX./(1-sX)).^(theta/(theta-1)) .*fx + k2 * ((sX>0).*(sD==1)) * FX + k2 * ((sX>0).*(sD<1)) * FXM ); % \varPsi
k7 = mean(phi.^(sg-1) .* (sX./(1-sX)) .*sD.^(-g*(sg-1)/(eps-1)));
k8 = mean(phi.^(sg-1) .* (1./(1-sX)) .*sD.^(1-g*(sg-1)/(eps-1))); 

P = (sg/ (sg-1))^(1/(1-g)) * (1-g)^(-1) * (g*beta^(eps/(eps-1))*qD)^(-g/(1-g)) * (mean(phi.^(sg-1).*sD.^(-g*(sg-1)/(eps-1))))^(1/((1-sg)*(1-g))); % This is P/w
bundle = (1-zeta)/zeta / (1-vartheta);
varrho  = 1 - (1/param.shareN-1)^bundle * (param.aN/(1-param.aN))^(vartheta*bundle) * (AN*P)^((zeta-1)/zeta) * (1- param.shareM)^(1/zeta);
w = pM / P * ( ((1-varrho)/varrho)^zeta * param.shareM/(1- param.shareM))^(1/(zeta-1));
sM = 1 / (1 + ((1-varrho)/varrho)^zeta * (P/(pM/w))^(1-zeta));
sN = 1 / (1 + ((1-param.aN)/param.aN)^vartheta * (AN * (varrho^zeta * (pM/w)^(1-zeta) + (1-varrho)^zeta * P^(1-zeta))^(1/(1-zeta))  )^(1-vartheta) );
aT = (1 - sN) * (1-sM);
Timplied = (  ( (1-sN)/aT * (1 - (sg-1)/sg * g *k5 *k8 * P^((sg-1)*(1-g))) - (1 + (sg-1)*(1-g))*(1/sg + k5 * k7 * P^((sg-1)*(1-g)) /sg) )/( (1 + (sg-1)*(1-g))*(1/sg + k5 * k7 * P^((sg-1)*(1-g)) /sg) + sN/aT * (1 - (sg-1)/sg * g *k5 *k8 * P^((sg-1)*(1-g))) ) )*param.L;
end